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ABSTRACT 

We present a phcnomcnological model of the dynamics of buoyant bubbles in the at- 
mosphere of a cluster of galaxies. The derived equations describe velocity, size, mass, 
temperature and density of the buoyant bubbles as functions of time based on sev- 
eral simple approximations. The constructed model is then used to interpret results 
of a numerical experiment of heating of the cluster core with buoyant bubbles in a 
hydrodynamical approximation (i.e. in the absence of magnetic fields, viscosity, and 
thermal diffusion). Based on the model parameters we discuss possible limitations of 
the numerical treatment of the problem, and highlight the main physical processes 
that govern the dynamics of bubbles in the intracluster medium. 

Key words: galaxies: cooling flows - galaxies: nuclei - galaxies: active - galaxies: 
clusters: general - galaxies: clusters: individual: Virgo - methods: numerical 



1 INTRODUCTION 

Depressions of the surface brightness in X-ray images of clus- 
ters of galaxies have been identified as low density and high 
temperatur e plasma bubbles (see, e.g., iMcNamar a et ail 
l200d , l200ll : iFabian et al.ll200rj| ), created by outflows from 
the central active galactic nucleus (AGN). Due to their high 
heat content these bubbles are thought to be one of the main 
heating sources of the cluster core, and the key ingredie nt 
in the solution of the cooling flow problem (|Fabianl 1994 ). 

The physical process of the deposition of the heat from 
the bubbles into the ambient ICM is not well understood. It 
is partly because we do not know all the relevant physical 
pro perties of the ICM (i t s turb ulent velocities, viscosity, etc., 
see ISchekochihin et al.l ([2005) for further discussion), and 
partly because theoretical models based on the numerical 
simulations with all the relevant physics are very complex, 
and more sophisticate d numerical models, e.g. including re- 
alistic magnetic fields l|Ruszkowski et aill2007l ). are not yet 
very common. 

In order to get a better understanding of the physi- 
cal processes that determine the behaviour of AGN-blown 
bubbles in the ICM we start from a basic hydrodynamical 
model. In the present article we analyse the physics behind 
the results of numerical simulations of the evolution of buoy- 
ant bubbles in the atmosphere of a cluster of galaxies in the 
absence of magnetic fields, viscosity, and thermal diffusion. 
We derive equations that describe velocity, size, mass, tem- 
perature and density of the buoyant bubbles as functions 
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of time based on several simple approximations. This ap- 
proach highlights the important hydrodynamical effects, and 
helps to understand and discuss limitations of the numerical 
framework for the description of AGN bubbles. It should be 
viewed as a first step in the modelling process, that can iso- 
late the phenomena caused by the simplest hydrodynamical 
effects. 

The structure of the article is as follows: in section [2] 
we outline the setup of the numerical simulation, which we 
later use to fit our model parameters; in section we de- 
scribe the vortex-ring model for AGN bubbles, derive the 
equations that determine the evolution of the parameters of 
the bubbles, and determine values for the free phenomeno- 
logical parameters of the model from comparison with the 
numerical data; in section [4] we analyse the physical prop- 
erties of the model, and how numerical artifacts can affect 
them; in section [5] we summarise the main findings and con- 
clusions. Appendices |A"1 [Bl and [Cl provide some background 
information about the hydrodynamics of the vortex rings. 



2 NUMERICAL EXPERIMENT 



In our previous work (|Pavlovski et al ] |2007l , hereafter PI) 
we have demonstrated that the evolution of buoyant bub- 
bles in the ICM is influenced by the hydrodynamical Kutta- 
Zhukovsky forces, which exist due to the circulation of the 
plasma in the space occupied by the bubble. One of the goals 
of the present work is to provide a quantitative description 
of the role of these forces in the overall dynamics of the 
bubbles. 
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Findings in PI were based on the analysis of numerical 
experiments of the buoyant ascent of bubbles in the ICM. 
The cluster atmosphere was m odelled using the ob serva- 
tional data for the Virgo cluster (|Ghizzardi et al ] |2004h . The 
bubbles were introduced ad hoc into an established cooling 
flow, at the point when the central temperature dropped be- 
low 1.5 keV. The initial temperature of the bubbles was fixed 
at T — 5 x 10 10 K (a factor of 10 3 larger than the tempera- 
ture of the ambient ICM). The density of the plasma inside 
the bubbles was calculated from the condition of the overall 
pressure equilibrium with the ambient medium (we tested 
two different pressure profiles). The initial bubbles' centres 
were placed symmetrically at the distance of 1.5ao from the 
centre of the cluster, where ao = 10 22 cm ~ 3.24 kpc was 
fixed as their initial radii, see Fig. [1] The size of the com- 
putational domain with periodic boundary conditions was 
10 24 cm. The computation was performed using the FLASH 
AMR hydrocode with 8 levels of refinement, resulting in a 
smallest cell size of w 4.88 x 10 20 cm. 

Our motivation for the choice of oo was the following. If 
E is the amount of energy released by an AGN, and Eq — fE 
is the fraction of the energy used to heat the ICM during 
the inflation of the bubble, then the volume of such a bubble 
is, 



Vo = E (i-l)/P, 



(1) 



where P is the ICM pressure, and 7 = 5/3 is the ratio of the 
specific heats. For typical AGN values this gives the size of 
the bubble as, 
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where L is an average AGN bolometric luminosity, and r 
is an average time between bubble injections, n am b is the 
number density of the ICM, T am b its temperature, and / 
sets the efficiency of the thermal coupling of the outflow 
with the ICM. We here assume / to lie in the range of a few 
percent (|Siiacki fc Springe]||2006h . 



3 PHENOMENOLOGICAL MODEL 

Analysis of the dynamics of large scale thermals and plumes 
in the atmosphere of the Earth (both nat ural and cre- 
ated during (nuclear) test explosions, see, e.g . , lTurner|[l957l ; 
lMorton|[l96a ; IWooddll997l : IZhou et al.ll200lf ) suggests that 
the evolution of hot bub bles in a stratified atmosph e re can 
be split into two stages. iHristianovich fc Rodionovl (|l954f ) 
proposed that during the first, relatively short stage of 
the ascent of a hot bubble it can be treated as remaining 
roughly spherical, while during the second stage its shape 
is best approximated by a torus. During the second stage 
the dynamics of the bubble is influenced by the velocity 
field around the bubble, which corresponds to that of a 
vortex ring. It is due to this velocity field that the mor- 
phology of the bubble changes (see PI for more details). 
Qualitative analysis of numerical simulatio of the hot bub- 
bles in the ICM in the absence of magnetic fields done by 
other research groups jGardinill2007l ; iRevnolds et"aT1l2005l ; 
iBriiggen fc Kaiserll2002 ) also support this picture. 




Xo=l-5a„ 



Figure 1. Initial size and position of the spherical bubble, and 
geometry and size of the torus-shaped bubble during the second 
stage of the ascent (not to scale). Point O marks the centre of 
the cluster, ao is the initial size of the bubble, x$ = 1.5ao is its 
initial position. R is the radius of the torus, r is the radius of 
the torus cross section, P is a point on the central circle axis line 
of the torus, vectors b and n are accordingly the bi-normal and 
the normal vectors of the central circle axis line at point P (see 
Appendix let . 



The morphological similarity of these (vastly different 
in size) bubbles is hardly surprising as the physical pro- 
cesses that govern them are very similar in the framework 
of a purely hydrodynamical model of the ICM and AGN 
bubbles. Such model is of course an oversimplification of the 
real phenomenon. AGN bubbles are filled with a relativistic 
magnetised plasma, which has an effect on the heat capac- 
ity of the plasma, and affects the diffusion properties at th e 
bubble/ICM interface (|Ruszkowski et al.ll2007l ; I Jonesll2007l ) . 
Here, however, we will assume that the magnetic fields are 
dynamically not important, and will view the plasma inside 
bubbles as thermal, simply having a much larger tempera- 
ture. 

In this rough approximation, both jet-blown ICM bub- 
bles and buoyant thermals in the Earth's atmosphere are 
characterised by high temperature and density contrast, be- 
ing, on average, in overall pressure equilibrium with the am- 
bient medium. Any departure from pressure equilibrium is 
smoothed out on the time scale of the order of ~ a/c, where 
a is the size of the bubble, and c = \J 'yP/p is the adiabatic 
speed of sound of the gas inside the bubble. The dynamical 
time scale, on the other hand, is ~ a/c am b, where c am b is the 
adiabatic speed of sound of the ambient medium. In the case 
of hot, low density bubbles, c am b/c <C 1, and we will assume 
that pressure equilibrium is maintained at all times. Finally, 
both hot bubbles in the ICM and hot atmospheric bubbles 
are rising due to the presence of the external gravitational 
field. 

Based on this qualitative comparison we present a phe- 
nomenological model of the evolut ion of the buoya nt bubbles 
in the ICM bas ed on the work of lOnufrievI (|l967t ) (for a re- 
cent review see iBelotserkovskv et al.l 2000l . in Russian) . 
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3.1 Growth of the spherical bubble 

In this section we will estimate the change of the size of a 
spherical bubble as it ascents in a stratified atmosphere of 
a cluster. In a static atmosphere the only physical processes 
that affect the volume of the bubble are the entrainment of 
the colder ambient plasma (at the same pressure), change 
of the ambient pressure (as the bubble ascents the ambient 
pressure falls) , and energy loss due to bremsstrahlung radia- 
tion. We can estimate the change of the volume due to these 
processes from the following simple argument. 

Consider a hot plasma bubble ascending buoyantly a 
distance dx during the time period dt, while a mass dM of 
the ambient gas at the temperature T am b is entrained into it. 
Note, that the entrained plasma needs to be accelerated to 
the velocity dx/dt = x to travel with the bubble. However, 
the fraction of the total energy of the bubble which goes into 
the increase of the kinetic energy of the entrained plasma is 
going to be much smaller than the fraction of energy which 
goes into heating of it to the new higher temperature T, 
provided that the ascent speed of the bubble is smaller than 
the speed of sound in the ambient ICM, jij < c am b- Here we 
will ignore the kinetic energy of the entrained mass. 

It is obvious that, 



dV_ 
V 



dM 
~W 



dp 
P ' 



(3) 



To estimate the right hand side of Q we note that in pres- 
sure equilibrium dp/p — — dT/T, and the temperature of 
the mixture, T, could be found from the equation of inter- 
nal energy conservation, 



T amb dM + TM = (T + dT)(M + dM). 



(4) 



where we have ignored differences in the heat capacities of 
the ambient plasma and plasma inside the bubble. Hence, 
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(5) 



where dpi is the change of density of the bubble due to 
the entrainment of the mass dM of the ambient plasma at 
temperature T am b, and £ is the contrast parameter, 



c = 



Pamb 



T 



(6) 



which can be related to the X-ray detectability of the 
bubbles; e.g., the ratio of the X-ray fluxes due to 
bremsstrahlung, oc p 2 T 1 ^ 2 , is equal to C 3 ^ 2 . 

Using equations © and © we get the following equa- 
tion for the change of volume due to the entrainment in the 
adiabatic homogeneous ICM, 



dV _ dM 
V ~ M 



(7) 



which has the simple solution, 
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(8) 



where Vo is the initial volume of the bubble, ao = 10 22 cm = 
3.24 kpc its initial radius and Co is the initial contrast pa- 
ramet er. Equation (JSj is due to iHristianovich fc Rodionovl 
(1954) . 

In our simulations bubbles are injected with the To ss 
10 10 K, and Co ~ 10 -4 . When the temperature of the bub- 
bles falls to, e.g., T = 2 x 10 7 K, C = 0.5, the corresponding 
relative change of the size of the spherical bubble (due to 
the entrainment) is a/ao « 1.26. It is worth noting that we 
have implicitly assumed here that the entrained material is 
mixed fully on time scales shorter than |ao/i|, and the bub- 
ble is homogeneous at all times. This assumption is generally 
not correct. The boundary layer of the bubble is entraining 
material, which is then getting mixed via turbulent motions 
or diffusion with the plasma in the bubble's centre. In this 
case the temperature of the bubble in not homogeneous: the 
temperature at the boundary of the bubble is lower then the 
temperature at the centre. When determining the size of the 
bubble from observations (or numerical data) , depending on 
an ad hoc threshold in some gas property used to find the 
boundary of the bubble, it can appear that the bubble is get- 
ting smaller, rather then larger, as ambient ICM is entrained 
into it. 

It is straightforward to show that the change of size of 
the bubble due to the change of the ambient pressure, and 
the loss of energy by bremsstrahlung radiation is small - 
generally not exceeding ten per cent. 

As the bubble ascents a distance dx, it reaches a layer 
with lower ambient pressure and adiabatically expands, 



dp 2 =p— . 



(9) 



Assuming pressure equilibrium and using equation Q we 

get, 



V_ 

Vo 



Po 



(10) 



For a distance x = 3.25ao travelled by the bubbles in our 
numerical experiment (see PI), we get a/ao ~ 1.12. 

The plasma inside the bubble emits bremsstrahlung ra- 
diation and cools. The density inside the bubble changes 
dp/dp = —dT/T = —dE/E, where dE is the loss of internal 
energy due to bremsstrahlung. The resulting change of the 
volume is, 



V_ 

Vo 



= 1 - 



AE 
~E~' 



(11) 



where AE is the total loss of internal energy due to 
bremsstrahlung during the ascent of the bubble. We can 
estimate the upper limit of AE/E using the values of den- 
sity and temperature of the bubble near the end of its as- 
cent. Using data from the numerical experiment (see PI): 
At « 7.7 x 10 14 s, p « 10" 26 gem -3 , T « 10 8 K, we get 
0.99 < a/a < 1. 

All of the mechanisms described above change the size 
of the bubble isotropically, and as a result the radius of the 
bubble towards the end of the ascent a « 1.35ao- To ac- 
count for the preferential sideways enlargement of the bub- 
bles, eventually leading to formation of a torus, and reconcile 
this analysis with the numerical data, we need to consider 
the main forces acting on the bubble during the ascent. 
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Figure 2. Right panel: forces acting on a segment (element) of the vortex ring. F b — buoyancy, F^ - drag, F^ z - Kutta-Zhukovsky force, 
_F kz ||, Fk z ± — components of the Kutta-Zhukovsky force. Left panel: velocity components. U - total velocity vector of the centre of mass 
of the vortex ring element, n - vertical component of the velocity vector, u± - horizontal velocity (expansion) due to the component 
Kutta-Zhukovsky force Fj sx ±. 



3.2 Buoyant Vortex Ring 

Sectio ns 13.21 13.31 13.41 13.61 are based on the w ork of lOnufrievI 
(| 19671 ') and lHristianovich fc Rodionovl l|l954h . 

The main force that determines the ascent of the AGN 
blown bubbles in the ICM is created by buoyancy due to 
the difference in densities of plasma inside the bubbles and 
the ambient ICM. During the ascent the mass entrainment 
creates friction forces which lead to the formation of a cir- 
culational motion of gas inside the bubble and around its 
surface. As the speed of the ascent grows, the spherical sym- 
metry breaks down. The initially spherical bubble expands 
sideways, since the pressure on the "top" and the "bottom" 
of the bubble becomes larger than the pressure on the sides 
(the velocity of the flow is larger on the sides). This process 
leads to the change of the shape of the bubble from spheri- 
cal into torus- like, see Fig[T] Inside the torus the gas rotates 
around the central circle axis line, and outside the torus a 
circulation flow is formed. 

The circulation of plasma around an element of the 
torus results in a Kutta-Zhukovsky force, Fk z , see Fig. [2] 
which is perpendicular to the overall velocity, U, of the el- 
ement. One of the components the Kutta-Zhukovsky force, 
Fk z x, is perpendicular to the direction of the ascent (_L- 
direction), and acts to enlarge the torus radius R. Another 
component of the Kutta-Zhukovsky force, F kz ||, is pointing 
in the direction opposite to the direction of the ascent (||- 
direction) and reduces the velocity of the ascent. 

The drag due to the viscosity and the kinetic energy 
of the wake underneath the bubble results in the drag force 
Fd, which also acts to slow down the ascent. Note taht the 
Fa acts in the direction opposite U, and hence not directly 
opposite to Ft,. 



3.3 Equation of motion 

The rate of change of momentum of an element of the vor- 
tex ring plus the associated rate of change of momentum 
of the surrounding medium (the latter is modelled as the 
added-mass, see Appendix [B| is equal to the sum of the 



forces acting on the vortex ring element, see Fig. f2] We can 
approximate the drag force using the formula applicable for 
a cylinder with a circular cro ss section of radius r moving 
through a fluid (|TrittorJll98Sl '). 

F d = ^C dPamb U\r 2 , (12) 

where Cd is a constant, p am b is the ambient density, U is 
the velocity of the vortex ring element, and r is the radius 
of its cross section. The experimental value of the constant 
Cd varies in range 0.1 .. . 100, and depends on the Reynolds 
number, Re, of the flow. Generally it is larger for flows with 
smaller Re. It is difficult to predict its value for our case, 
since the "cylinder" in this case is an element of the vortex 
ring, which is not a solid body, as it is under-dense com- 
pared to the ambient medium. The main motivation for us- 
ing equation (|12[l is purely geometrical. 

The Kutta-Zhukovsky force is proportional to the cir- 
culation, r, (see appendix [Cl and PI) of the flow outside the 
vortex ring, and is perpendicular to the velocity, U, of the 
element of the torus (vortex ring), 

Fkz = Tp amb URda, (13) 

where R is radius of the torus, and a is the central angle. 

By projecting the vector quantities on the || and _L di- 
rections, we get the following equations of motion, 

^ (tt|| (M + M')) = gV Pamh (l - <) - 27rp amb rFu ± 

1 2 ( 14 ) 
— 2 < ^ d7rr 

and 

^ («i (M + M')) = -2 7 rp amb riiM|| - i(7d7rrVt/, (15) 

where M' = 2n 2 r 2 Rp am h is the added-mass for the torus 
(see Appendix [B}, g is the gravitational acceleration, and 

U = ^Ju\ + (16) 

is the velocity of the torus element. 
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3.4 Heat balance 



3.7 Circulation 



Entrainment of the cooler ambient plasma changes the den- 
sity of the bubble according to equation ((5J , while the bub- 
ble also expands due to the fall of the ambient pressure as 
given by equation (|lip . From the definition of the contrast 
parameter, 

d£ _ dp _ d/Oamb 
C P Pamb 

we get the following formula for the change of contrast due 
to the mass entrainment and the change in the ambient pres- 
sure and density, 

dpamb 



dC ,,dM dP 



Pamb 



(17) 



which implies the following differential equation, 



dt 



= c 



(i-C) 



1 AM 

M~dT 



1 dP 

jP dx 



1 dp am b 



Pamb 



dx 



(18) 

The last term in (I18|l . which is enclosed in round paren- 
thesis, is a measure of the deviation of the ambient atmo- 
sphere from adiabatic conditions. For an adiabatic atmo- 
sphere dp am b/p a mb = dP/ (-P7), and the term is exactly zero. 

3.5 Mass entrainment 

The entrainment of the ambient ICM into the bubble hap- 
pens due to growing fluid instabilities on the surface of the 
bubble. The boundary layer of the buoyant bubble is corru- 
gated by the Rayleigh- Taylor (RT) instability and becomes 
turbulent. The interpenetration distance h of the light fluid 
into the heavy fluid in a classical RT instability scenario can 
be written as (see, e.o.. lShardfl984l ). 

h = C^Atgt 2 , (19) 

where M = (p amb - p)/(p am b + p) = (1 - C)/(l + C) is the 
Atwood number, g is gravi tational acceler ation, t is time, 
and Ct is a constant (see also lLiu et al .120071 , for a discussion 
of numerical methods for mass diffusion and values of the 
constant). Using (|19[) we get the following equation for the 
mass transfer rate due to the RT instability, 



1 dM 1_ 

M dt m l + (r 



(20) 



were we have substituted U for gt as it is a more appropriate 
measure of the velocity of the bubble/ICM interface in our 
case. 

3.6 Size of the bubble 

Using the formula for the mass of the torus, 
M 



2n 2 Rr 2 p 



(21) 

it is straightforward to write an equation for the change of 
the sizes, R and r, of the torus, 



dr r 
dt = 2 
or, equivalently, 
d_R 
dt 



1 dM 
M~dT 



= u±, 



dr 
dt 



1 dM 

M~dT 



i_dR 
R~dt 



u± 
R 



C dt 



C dt 



dp a 



Pamb da; 



dp a 



Pamb 



dx 



(22) 



(23) 



The circulation of the flow around the element of the torus is 
determined by the angular velocity of the fluid inside it. As 
the mass of the torus grows, the angular velocity of the fluid 
inside it decreases. Any disconnection between the inner and 
the outer circulation flows (i.e. difference in their angular ve- 
locity) leads to the production of vortices near the surface 
of the torus, which then are shredded away, slowing down 
the rotation of the flow further. It is reasonable to assume 
that the angular velocity is directly proportional to the cir- 
culation of the flow around the element of the torus, uj oc V. 
Consider a cylinder of radius r rotating in a fluid with kine- 
matic viscosity v aroun d its axis of symmetry. T he torque 
due to viscous stresses (lLandau fc Lifshitz| [l987. §18) per 
unit length of this cylinder is 87r 2 i/p am b^?r 2 ciJ = AvMui/Q. 
The torque of the friction forces will change the angular 
momentum of the cylinder, l/2Mr 2 u. This implies that, 



— (MrT) oc vMT/C- 
dt 



(24) 



In the case of the vortex ring the role of the viscosity (friction 
forces) is played by the interpenetration of the fluids across 
the boundary of the torus. The movement of the plasma 
across the boundary, in other words the rate of mass en- 
trainment, determines whether the motion of the fluid on 
one side of the boundary has an affect on the motion of 
the fluid on the other side of the boundary. Formally, from 
dimensional arguments we can set u oc (^r 2 M~ 1 dM/dt. Sub- 
stituting the expression for v into equation (|24p and assum- 
ing that the rate of change of the radius, r~ 1 dr/dt, is much 
smaller than the rate of change of the mass, M _1 dM/di, we 
get the following equation for the change of circulation, 



1 dr _ 1 dM 
T~dt ~ ~ c M"dT' 



(25) 



where C c is a constant which needs to be determined from 
comparison with the data. 



3.8 Simultaneous ODEs 

The velocity of accent x is the sum of the vertical velocity of 
the element, tt||, which is determined by the buoyancy and 
the hydrodynamical forces, plus the self-induced vortex ring 
velocity, Wind) which is the result of the curvature (1/R) of 
the vortex ring (see Appendix [C] for more information), 



dx 

- = u ll+ v ind , 



(26) 



where the induction velocity «i n d is proportional to the cir- 
culation and the curvature of the vortex ring, oc T/R. In the 
present model we will use the expression for Wi n< j given by 
equation (|C 10|) in Appendix [C] 

Together equations ([26]), (J25J) , ([23]), (3TJI, {20]), (IB). 
(|16|) . (|15p . (|14l) form a closed system of ordinary differen- 
tial equations (ODEs), which determines the values of eight 
physical parameters, u±, uu, R, r, x, F, M, and £ as func- 
tions of time. Three phenomenological parameters, Cd, C m , 
and C c , need to be determined by comparison with the nu- 
merical data. 
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3.9 Solution 

The simultaneous ODEs in their dimensionless form are 
written in full in Appendix [X] To determine the values of 
the phenomenological constants, C C; d,m, we have fitted the 
solution of the ODEs to the data from our numerical exper- 
iment C described in PI. As the initial condition we have 
selected a point when the bubble shape changes from spheri- 
cal to toroidal. In the simulation this happens approximately 
16 Myr after the injection of the bubble, and to = 722 Myr 
after the start of the simulation. For reference, the dimen- 
sionless time unit is approximately 8.1 Myr, the velocity unit 



is 3.9 x 10 7 kms _i , the distance unit is 10^ 2 cm, and the 
mass unit is 9.8 x 10 40 g. The initial values of the physical 
parameters at to were measured using the averaging tech- 
niques described in PI. Their numerical values were found 
to be as follows, u±(to) = 0.26, U||(to) = 0.42, R(t Q ) = 0.99, 
r(t ) = 0.72, x(t ) = 2.56, r(t ) = 2.71, M(t ) = 2.4x 10" 2 , 
C(*o) = 2.4 x 10" 3 . 

We solve the ODEs numerically using the adaptive step 
Bulirsch-Stoer method (|Press et alj|2002h . To find values of 
the phenomenological constants we minimised the function, 



22 




Figure 3. Growth of the bubble mass: solid line is the best fit 
model, points with vertical error bars are measurements from the 
numerical simulation. 



EE 

3=1 i=l 



(27) 



where tj are times of the snapshots of the computational 
grid produced during the numerical simulation (we have 42 
snapshots in total), A; = {M, R, r, x, f } are the parame- 
ters of the bubble (torus) calculated from these snapshots, 
AXi are their estimated errors, and Yi(tj) are the param- 
eters of the bubble (torus) from the solution of the simul- 
taneous ODEs for the times tj. The estimated errors, AX;, 
were calculated from one sigma variations of the statisti- 
cally determined quantities (see PI), and errors due to the 
finite grid resolution. We did not include into Xi parame- 
ters that rely on the estimation of velocities, since the com- 
paratively large and varying time between the snapshots, 
t(J) = tj — tj-i 7^ const, makes it difficult to determine the 
values of the velocities accurately. 

The best fit values for the constants were found to be 
C c = 0.01, Cm = 4, Cd = 24. Due to the uncertainties 
in the values of the parameters at to (especially velocities 
and circulation) the resulting values of the constants are not 
determined precisely. The values of the constants are going 
to be different if different initial conditions are used. By 
varying the initial parameters by ten percent, we found that 
the corresponding variation of the values of the constants 
was approximately limited to the ranges, C c ± 0.05, C m ± 1, 
C d ± 10. 



4 ANALYSIS OF SOLUTION 



In Figs. [3] |4][5l and[6]we compare the evolution of the param- 
eters determined from the numerical simulation, Xi, with 
the prediction of our model, Yi. It is clear that the model 
with the best fit constants determined above, reproduces the 
data from the numerical experiment quite well. We achieve 
a good fit for the contrast parameter and the mass of the 
bubble over three orders of magnitude in scale. The fit for 
the bubble radii seems to be a bit more problematic. It is 



1 .000 F 



0.100 r 



0.010 r 



0.001 




Figure 4. Change of the contrast parameter with time: solid line 
is the best fit model, points with error bars are measurements of 
the contrast paranter from the numerical simulation. 



important to note, however, that by changing the initial con- 
ditions for the solution of the ODEs we can get rather better 
fits for radii (at least by eye). This sensitivity to the initial 
conditions is not surprising. We find that all Yi, but espe- 
cially radii and position, depend on the value of the initial 
circulation, T(to). Since the constant C c « 0, the circula- 
tion remains roughly constant in time (conserved), and the 
precision of the measurement of the initial value is reflected 
in the overall quality of the fit of the parameters that de- 
pend on it directly, i.e. R, r, and x; the latter through the 
self-induced velocity of the vortex, v,^. 

To illustrate the matching of the model to the three- 
dimensional distribution of the bubble material at different 
times as computed in our simulation we compare them di- 
rectly in Figs. [7] [8] and[9] Note that the approximation of the 
bubble shape with a torus at to, Fig. [7] is not perfect. The 
distribution of the bubble material, although resembling a 
torus, does not have a perfectly circular cross section. The 
cross section of the bubble's torus is elliptical, and the size 
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Figure 5. Growth of the bubble: solid lines are the best fit mod- 
els (upper curve for R, lower curve (red) for r), dashes with the 
vertical error bars are measurement from the numerical simula- 
tion. 
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Figure 6. Ascent of the bubble: solid line — the best fit model, 
dashes with error bars - measurement from the numerical simu- 
lation. 



of the opening at the top of the bubble is smaller than the 
size of the opening at the bottom of the bubble. At later 
times the matching of the shape of the bubble to a torus 
with a circular cross section is much better. Despite that, 
our model was able to capture the overall size and shape of 
the bubbles quite well. 

The solution of the ODEs suggests that the vertical ve- 
locity of the bubbles due to the hydrodynamical forces, un, 
quickly falls due to significant drag, as follows from the high 
value of Cd ~ 30. In fact, we find that the ascent veloc- 
ity of the bubble, x, in the current model is dominated by 
the self-induced velocity of the vortex ring, wind, which is 
an order of magnitude larger than the velocity mm. A large 
drag coefficient is generally a characteristic of low Reynolds 
number (high viscosity) flows. Although we did not model 
any physical viscosity in our numerical simulation, the resid- 
ual numerical viscosity is always going to be present. It is 
difficult to quantify its contribution. Qualitatively, however, 



we can state that the velocity of the ascent is more signif- 
icantly affected by the the self-induced velocity of the vor- 
tex in the medium with high viscosity, and therefore large 
drag. In the framework of this model the bubbles will never 
reach "a ceiling" i.e. a maximum height. They will deposit 
their thermal content along their ascent path, but even af- 
ter they reach £ = 1 and disappear the residual vortex will 
keep moving towards the outskirts of the cluster 1 . If the 
ICM has very low viscosity, the Reynolds number is likely 
to be much larger than a value of few hundred, which is 
characteristic for our present simulation. In this case the co- 
efficient Cd can be much smaller, and x will no longer be 
determined by the self-induced velocity, but rather by the 
balance between the buoyancy and the vertical component 
of the Kutta-Zhukovsky force. 

Another possible pitfall for our model is the overesti- 
mation of the mass diffusion across the surface of the bub- 
ble from the fit to the simulation data. It is well known 
|Liu et all l2007h that generic hydro schemes tend to be 
super-diffusive, and yield too large values of coefficients like 
C m . We can roughly estimate the impact of a possible reduc- 
tion of the coefficient in the following way. By substituting 
equation (|20[) into equation (|18|) . and neglecting the proper- 
ties of the atmosphere we get, 



dC 1 + C 

<(1-C) ; 



= a 



u 



(28) 



which is easily integrable if we assume that U and r are 
constants, 



2(C - Co) 

(i-C)(i -Co) 



(29) 



The left hand side (l.h.s.) of equation (|29[) is a finite number 
for any given values of the contrast parameter < £o < C < 
1 (l.h.s. w 27 in our case: £ = 0.9, Co = 0.001). Therefore, if 
the value of C m is reduced by a certain factor, /, it will take / 
times longer for bubbles to reach the same contrast f. Given 
that the typical mass diffusion constant in simu lations of 
RT in stability is overestimated by a factor of two (|Liu et al.l 
l2007h . it is quite likely that the bubbles in the ICM survive 
at least twice as long as a typical simulation would tend to 
suggest. 

The small value of the constant C c w implies that the 
circulation in our model is approximately conserved, V ~ 
const, and the fluid is nearly ideal and incompressible. The 
opposite extreme C c = 1 corresponds to the conservation 
of angular momentum, FM — const, see (I25|l . and would 
describe a bubble unaffected by the surrounding medium. 



5 DISCUSSION AND SUMMARY 

The developed model is by no means a full description of 
the physics of buoyant bubbles in real clusters of galaxies. 
In this work we have tried to investigate the simplest pos- 
sible approximation, and understand how different physical 



1 The drag coefficient in problem of the rise o f the cloud from 
a nuclear explosion is much smaller, Cd ~ 0.5. Onufrieyj dl967T ) 
found that in this case the cloud reaches the maximum height, 
which depends on the stratification of the atmosphere and its 
initial size. 
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Figure 8. Comparison of the distribution of the bubble material and the approximating torus at t = 93.9. 



properties of the bubbles are interconnected. Hydrodynam- 
ical approximations have been widely used to approximate 
processes of heating in clusters with AGN bubbles (e.g., 
iBriiggen fc Kaiseill2002l ; lGardinJl2007r i, and it is important 
to understand the properties of such solutions in order to 
know their limitations. Magnetic fields are clearly impor- 
tant ingredients in AGN physics. They regulate dep osition 
of heat from the bubble into the ambient medium (|jonesl 
120071 ; iRuszkowski et al]|2007ft by suppressing thermal con- 
duction around the surface of the bubble, which otherwise 
can lead to a rapid evaporation of the bubble (see PI). 

To summarise, in this study we have developed a phe- 
nomenological model, which describes the dynamics of buoy- 
ant AGN bubbles in the atmospheres of clusters of galaxies. 
The value of the free phenomenological constants was de- 
termined from a comparison with the simulation data. The 
resulting fit of the prediction of the model to the simula- 
tion data is good. We have analysed possible outcomes that 



correspond to different values of the constants. The main 
points are as follows: 1. In low Reynolds number flows the 
drag coefficient is large, and the overall dynamics is deter- 
mined by the circulation of the flow around the bubble. 2. 
Since mass entrainment is likely to be overestimated in nu- 
merical simulations, the bubbles can potentially survive in 
the ICM much longer than predicted by a typical simula- 
tion. 3. The conservation of circulation in our model is a 
consequence of the near incompressibility of the ICM (given 
the characteristic velocities), and absence of viscosity. Note, 
that the numerical viscosity alone may significantly decrease 
the flows Reynolds number. 

The model presented here provides a framework for the 
interpretation of numerical and observational results. This 
purely hydrodynamical basis can now be extended to include 
more of the important physics, e.g. the effects of magnetic 
fields. 
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Figure 9. Comparison of the distribution of the bubble material and the approximating torus at t = 99.4. 
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APPENDIX A: SYSTEM OF ODES 

Equations (25), (25), (23), (2T), (20), IQD, (3D. G3>> and 
(|14p , which were derived in sections § 13.21 - 13.81 form a sys- 
tem of ordinary differential equations (simultaneous ODEs). 
The solution of these equations determines values of the bub- 
ble sideways enlargement velocity u±, vertical velocity due 
to hydrodynamical forces un, radial size R, vertical size r, 
position x, circulation F, mass M, and contrast parame- 
ter £ as functions of time. The identified physical processes 
were parametrised with the phenomenological parameters 
for drag Cd, mass entrainment C m , and circulation C c . 

It is convenient to work with dimensionless variables, 
We since measure sizes and distances in terms of the initial 
size of the bubble ao . We can use factors ^/aogo for velocities 
and s/ao/ go for time to make these variables dimensionless, 
where go is the gravitational acceleration at the starting 
point go = g(xo). By substituting, 



{r,R,x} 



{r,R,x} 



ao 
t 



{x,u±,uu} 



{x,u ± ,u ll } 
y/aogo 



ao^/aogo 



(Al) 



Ar _ r f" 1 AM u± 1 AC, 
dt ~ 2 [M~dl ~R ~ C dT ~ 
dR 

dT =U± ' 
its characteristic velocities, 

du|| _ 1 — C Tu± 
~~aT = 1 + ( 9 ~ 7r(l + C)r 2 

1 dC 1 dM 



x— log p amb 
ax 



(A3) 



+ «n 



C d 



u 



C(l + C) dt M At 4tt R(l + 



Tu 



Au± 

~~aT ~ tt(1 + C)r 2 



(A4) 



+ Mj 



1 



1 dM C" d 



U 



C(l + dt M dt 4tt + C) 



its position (height), 



a; = U|| + Wind, 

r 



Wind = 



4ttR 



, 8R 1 
log - 

r 4 



(A5) 



121og(8i?/r) - 15 / r - - 
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circulation of the ICM around it, 



(s) 



dr _ r dM 

dt ~~ ~ c M~dT' 



and the mass entrainment rate, 



AM 
At 



CUM 



1 + C r 



(A6) 



(A7) 



into the original equations, we get, after some simple al- 
gebra, the following ODEs describing the evolution of the 
parameters of the bubble: its contrast, 



dt s 



^-^m^t +x a-x 1o ^ p/ ^ 



(A2) 



APPENDIX B: ADDED-MASS 

A body moving through a continuous medium creates a ve- 
locity field in that medium. A spherical bubble ascending in 
a static atmoshpere of a cluster of galaxies creates a velocity 
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field in the cluster. The total kinetic energy of the system, 
which consists of the moving bubble plus the ambient at- 
mosphere, is the sum of the kinetic energy of the of the 
medium (the kinetic energy of the induced velocity field), 
plus the kinetic energy of the bubble. It is possible to show 
that such motion can be described as a motion of a spheri- 
cal bubble with an effective mass, which consists of its own 
mass plus an added-mass, which is proportional to the mass 
of the medium displaced by the body. Here we will assume 
that the flow is laminar, potential (i.e. non-rotational), and 
the bubble has a constant radius. 

The velocity field in an ideal incompressible fluid with 
zero vorticity can be described by the equation, 



V x v = 0, 



(Bl) 



or equivalently, 

v = -Vy>. (B2) 
Using the incompressibility (V ■ v = 0) condition we get, 



V> = 0, 



(B3) 



which is the well-known Laplace equation for the velocity 
potential, ip. 

Consider a uniform flow of an incompressible, non- 
rotational fluid past a sphere. By uniform flow here we mean 
that at large distances from the sphere the flow has a uni- 
form velocity C/k, where k is the unit vector (along the z 
axis) of a Cartesian coordinate system. We can chose the 
centre of the sphere as the centre for the coordinate sys- 
tem. If the radius of the sphere is a, then we have to solve 
the Laplace equation (|B3[) in the region a < r < oo. The 
boundary condition applied to the surface of the sphere is, 



dr 



= 0. 



(B4) 



The boundary condition at infinity can be taken to be, 

tp — —Ur cos 6, at r — > oo, (B5) 

which would give the velocity (in a spherical system of co- 
ordinates) , 

v(r, 0,<j>) — U cos 9e r — U sin 6eg, when r — > oo, (B6) 

where we used the relation k = cos 9e r — sin 8eg ■ 

The general solution of the Laplace equation ()B3|I in 
spherical coordinates is, 



oo I 



ip =J2 [Air 1 + B ir 



Pi m (cos6>)x 



1=0 m=0 

x [S™ sin <pm + C™ cos 4>m] , 



(B7) 



where p m are the associated Legendre functions 
(|Morse fe Feshbachlll953h . 

The boundary conditions (|B4|) and (|B5[) can be satisfied 
only for the following combination, 

tli = -u(r+^\cos6. (B8) 



The velocity of the flow can be easily found using ()B2[) . 

v = (7k — U^-r ( cos 9e r — - sin 6eg 
r 6 \ 2 



(B9) 



With a simple Galilean transformation we can now con- 
sider the problem of motion of the sphere with velocity t/k 
through a fluid at rest (at infinity). If we take the origin of 
the new coordinate system at the instantaneous position of 
the centre of the sphere, then the flow pattern around the 
sphere is given by, 

„3 



-U- 



cos 6e r 



1 



sin 8et 



(B10) 



The kinetic energy of the fluid around the sphere is, 



1 r 2 * r r°° , 2 2 

A"fiuid =-p / v r sin0drd0d<p 

2 Jo JO Ja 
1 3 rr 2 

= -Tvpa U 

=-M'U 2 , 
2 



(BH) 



where M' = 1/2 Vp, V = 4/37ra 3 . The total kinetic energy 
of the system fluid plus the sphere is equal to the sum of the 
kinetic energy of the sphere and the kinetic energy of the 
fluid, 



K = - (M + M') U 2 



(B12) 



and can be interpreted as the kinetic energy of the sphere 
with the effective mass of M + M' - the mass of the sphere 
plus an added-mass equal to half of the displaced mass of 
the fluid. 

A similar d erivation for a cylinder yields (see e.g., 
IChoudhurilll998l ) M' = Vp, where V is the volume of the 
cylinder. If we consider an element of the vortex ring to 
be roughly cylindrical, then for a whole torus moving in 
the || direction, the added-mass also should be equal to 
the total mass of the displaced fluid, M' = Vp, where 
V = 2-K 2 r 2 R is the volume of a torus. A rigorous mathe- 
matical treatment of the problem (|Miloh et al.lll978T ) proves 
this approximation to be accurate to within a few per cent, 
1 < M' /{Vp) < 1.0625. 

Note that the motion of two identical spheres in op- 
posite directions along a line connecting their centres with 
velocities U = Uj = —U i, yields the expression for the ki- 
netic energy of the fluid l|Lambll 19321 ). 



2K = pV 1 



3 a A 
16 



+ o 



(B13) 



where V = 4/3na 3 is the volume of the sphere, and x is half 
of the distance between the spheres. Equation (|B13[) implies 
that the resulting effective mass for a sphere in the presence 
of an identical sphere moving in the opposite direction is 
slightly lower than that of a single sphere. Therefore, our 
approximation of the added-mass of the torus, M' — Vp, 
which is slightly lower than the va l ue for a single torus as 
given by calculations of Miloh et alj (Il978l ) is, in fact, appro- 
priate since in our cluster there are two bubbles, i.e. there 
is a second torus, which is moving with the same speed in 
the opposite direction. 



APPENDIX C: VELOCITY INDUCED BY THE 
VORTEX RING 

The motion of the gas in cluster around the toroidal bubble 
is not irrotational. While the results of the previous section 
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remain valid, they are incomplete. To describe motion of the 
gas around a toroidal bubble we have to take into account 
the rotation of the gas around the centre line of the torus. 
The rotation itself does not change the value of the added- 
mass, but it is essential for proper treatment of the velocity 
field around the toroidal bubble. 

The following genera l discussion is based on the mate- 
rial from iBatchelorl (|l967T ) . 

Consider an incompressible fluid with the velocity field, 

v, 

V ■ v = 0, 

which can be described in terms of the vector potential, B, 
v = VxB, (CI) 

and vorticity vector ui, 

uj — V x v. 

The above equations imply 



V(V-B)-V 2 B = , 



(C2) 

Se can seek a solution of (|C2|I in the volume where VB = 0, 
so that, 



and 



B(x) 



4* J 



dU(x') 



(C3) 



where s = jx — x'|. Using (IC3[) the gauge V • B = reduces 
to w • n = 0, where n is the normal vector to the boundary 
of the fluid. The latter can always be satisfied in a volume 
possi bly extended beyond the physical boundary (jBatchelorl 
Il967l . §2.4). Using (JCU we find the velocity of the fluid in 
terms of the vorticity vector, 



v(x) 



1 

4% 



s x w(x') 



dU(x') 



(C4) 



The vortex-line (by analogy with the flow line) is the line 
whose tangent vector is everywhere parallel to the local value 
of the vorticity w. The vortex-tube is a surface constructed 
by all vortex lines passing through a given reducible closed 
curve C 2 . The integral of vorticity over an open surface A 
bounded by the same closed curve C, lying on a vortex tube 
and passing round it once is independent of the choice of the 
line and, therefore, the surface A. In other words, the flux of 
vorticity through the vorticity-tube, along the vortex-line, is 
conserved. The value of the surface integral over the surface 
A, i.e. the flux of vorticity through a cross section of the 
vortex-tube, is called the strength of the vortex tube, 



r = / w ■ ndA = 



Ax. 



The line integral is taken over a closed line, C, lying on the 
vortex-tube and passing around it once. It is called circula- 
tion. 

2 This definition of the vortex-tube is identical to the definition 
of the flux tubes in a magnetic field. Further analogies do also 
hold, e.g., | |C6|| is analogues to the Biot-Savart law, where the 
line vortex plays the role of a current, and velocity plays the role 
of the resulting magnetic field. 



The line vortex is a useful mathematical abstraction in 
cases when the vorticity is large in a limited volume (in the 
vicinity of a line) and negligible elsewhere. If 81 is a vector 
element of the line vortex which lies in the volume SV, by 
the definition we have, 



Jsv 



udv = r<ji, 



(C5) 



where F is the strength of the vortex tube constructed 
around the line vortex. Using (|C5|) and (|C4|I it is straight- 
forward to get the velocity distribution of the fluid, 



v(x) 



4% 



s x dl(x') 



(C6) 



where s is the vector connecting the point x in the flow with 
the point x' on the line vortex. 



CI Application to buoyant bubbles 

We have argued (PI) that the velocity field of plasma in- 
duced by the ascending bubbles is that of a vortex ring, i. e. 
a circular line vortex of radius R and strength T. Using (IC6|I 
it is easy to get an expression of the velocity of the plasma 
directly under the rising circular line vortex (vortex ring), 



v(x) = 



YR 2 



2(x 2 + R 2 ) 3 / 2 ' 



(C7) 



where x is the distance from the centre of the bubble to the 
point directly under it. This is the velocity with which the 
material will be lifted up from the centre of the cluster by 
the ascending bubble due to the circulation. 

The velocity field induced by the vortex ring has an- 
other interesting property. It is obvious from (|C6P that there 
is a singularity in the velocity field at the points on the vor- 
tex. Careful exami nation of the vel ocity field in the vicinity 
of the line vortex (Batchelor 196/]) shows that the velocity 
field near the vortex consists of the rotational motion around 
it plus a translational motion. If t, n, b are the tangent, the 
normal, and the binormal vectors of the line vortex at the 
point P, then the position of a point in the plane perpen- 
dicular to the line vortex at P can be written as, 

r = q 2 n + q 3 h. 



The velocity field at the distance r = \/ a 2 + a\ — > near 
the point P, see Fig.[TJ has the asymptotics (|Batchelorlll967l . 
§ 7.1), 



■ (b cos i 



n sin 0) + -£-b log -+0 (r°) , (C8) 



27rr v T ' AttR 

where <f> is a polar angle in the plane defined by the vectors 



n and b. The first term in (|C8|I represents the expected 
circulatory motion about the line vortex. The second term 
shows that there is another weaker singularity of the velocity 
distribution associated with the local curvature of the line 
vortex, this is the induced velocity, 

V ind = bVind, 

The velocity of the fluid in the vicinity of point P has a large 
velocity in the direction of the binormal, with a magnitude 
varying asymptotically as log 1/r. 
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iLambl l| 19321 ) (Ch. 7) gives a proof of the Kelvin formula 
for the self induced velocity of a circular vortex ring of a 
small cross section in a perfect fluid, 



Wind 



r 

4-kR 



lot 



811 
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(C9) 



In case when the vorticity is confined to the surface of the 
ring (s o called hol low vortex), the term —1/4 is replaced by 
-1/2 |Hicksil 18831 ). 

Equation (|C9[) was applied to buoyant vortex r ings by 
Tumerf <|l957t > (see also lMortonlll960l ; lTurner|[l969l; I Woods 
1997). The general problem was re-examined by iFraenkel 
( 1972), who proved the existence of steady vortex rings, and 
gave a generalisation of (IC9|l f or arbitrary dis tributions of 
vorticity as a function o f r, an d Saffmanj (|l970P l for the case 
of viscous vortex rings. lOnufrievi (|l967l ) used the following 
formula for the induced velocity for description of motion of 
the cloud from a nuclear explosion, 



Wind 



r 

121og(8ii/r) - 



8R 



32 



Ti) 



(CIO) 



REFERENCES 

Batchelor G. K., 1967, An introduction to fluid dynamics. 

Cambridge University Press 
Belotserkovsky O. M., Andruschenko V. A., Shevelev Y. D., 

2000, Dynamics of spacial and vortex flows in the inhomo- 

geneous atmosphere: A numerical experiment. Moscow: 

Yanus-K 

Briiggen M., Kaiser C. R., 2002, Nat., 418, 301 
Choudhuri A. R., 1998, The Physics of Fluids and Plasmas. 

Cambridge University Press 
Fabian A. C, 1994, ARA&A, 32, 277 

Fabian A. C, Sanders J. S., Taylor G. B., Allen S. W., 
Crawford C. S., Johnstone R. M., Iwasawa K., 2006, MN- 
RAS, 366, 417 

Fraenkel L. E., 1972, J. Fluid Mech., 51, 119 

Gardini A., 2007, A&A, 464, 143 

Ghizzardi S., Molendi S., Pizzolato F., De Grandi S., 2004, 
ApJ, 609, 638 

Hicks W. M., 1883, Royal Society of London Proceedings 
Series I, 35, 304 

Hristianovich S. A., Rodionov V. N., 1954, Technical re- 
port, On the rise of a cloud. Inst. Chem. Phys. Ac. Sci. 
USSR 

Jones T. W., 2007, ArXiv e-prints, 708 
Lamb H., 1932, Hydrodynamics. Cambridge University 
Press 

Landau L. D., Lifshitz E. M., 1987, Fluid Mechanics, sec- 
ond edn. Vol. 4, Pergamon Press 

Liu X., Li Y., Glimm J., Li X. L., 2007, Journal of Com- 
putational Physics, 222, 644 

McNamara B. R., Wise M., Nulsen P. E. J., David L. P., 
Sarazin C. L., Bautz M., Markevitch M., Vikhlinin A., 
Forman W. R., Jones C, Harris D. E., 2000, ApJ, 534, 
L135 

McNamara B. R., Wise M. W., Nulsen P. E. J., David L. P., 



Carilli C. L., Sarazin C. L., O'Dea C. P., Houck J., Don- 
ahue M., Baum S., Voit M., O'Connell R. W., Koekemoer 

A. , 2001, ApJ, 562, L149 

Miloh T., Waismann C, Weihs D., 1978, Jour, of Eng. 
Math., 12, 1 

Morse P. M., Feshbach H., 1953, Methods of Theoretical 

Physics, Part II.. New York: McGraw-Hill 
Morton B. R., 1960, J. Fluid Mech., 9, 107 
Onufriev A. T., 1967, Jour, of Appl. Mech. and Thech. 

Phys., 8, 3 

Pavlovski G., Kaiser C. R., Pope E. C. D., Fangohr H., 
2007, Morphology of flows and buoyant bubbles in the 
Virgo cluster, submitted, arXiv:0709.1790 

Press W. H., Teukolsky S. A., Vetterling W. A., Flannery 

B. P., 2002, Numerical Recipes in C++. The art of sci- 
entific computing, second edn. Vol. Volume 1 of Fortran 
Numerical Recipes, Cambridge University Press 

Reynolds C. S., McKernan B., Fabian A. C, Stone J. M., 
Vernaleo J. C, 2005, MNRAS, 357, 242 

Ruszkowski M., Enfilin T. A., Briiggen M., Heinz S., Pfrom- 
mer C, 2007, MNRAS, pp 400-+ 

Saffman P. C, 1970, Stud. Appl. Math., 49, 371 

Schekochihin A. A., Cowley S. C, Kulsrud R. M., Hammett 
G. W., Sharma P., 2005, in 'The Magnetized Plasma in 
Galaxy Evolution' Krakow, Poland, Sept. 27th - Oct. 1st, 
2004 Magnetised plasma turbulence in clusters of galaxies, 
p. 200 

Sharp D. H., 1984, Physica D Nonlinear Phenomena, 12, 3 

Sijacki D., Springel V., 2006, MNRAS, 366, 397 

Tritton D. J., 1988, Physical Fluid Dynamics, 2nd ed. edn. 

Oxford, England: Oxford University Press 
Turner J. S., 1957, Proc. Roy. Soc. A., 239, 61 
Turner J. S., 1969, Annu. Rev. Fluid Mech., 1, 29 
Woods A. W., 1997, J. Fluid. Mech., 345, 347 
Zhou X., Luo K. H., Williams J. J. R., 2001, Eur. J. Mech. 

B. - Fluids, 20, 233 



